### This script was written by Eric Deal, 06.22.2020
### Script should be run with python 3

#%%
import numpy as np, pandas as pd

# Grain density was measured by weighing a sample of 15-40 grains in an empty 10
# ml vial. The vial was then filled with water to the 10ml mark, and the mass of
# the vial was remeasured. The mass of the vial when empty and when filled with
# 10 ml of water was known. From this we could obtain the volume of the grain
# sample by calculating the difference between the vial before and after it was
# filled with water and dividing by the known density of water.  The average
# grain volume per sample was then obtained by dividing this volume by the
# number of grains in the sample. The procedure was repeated 3 times for each
# grain type. The largest source of error was our ability to fill the vial to
# precisely 10 ml. We used a high-precision graduated vial with a pipette and
# estimate our ability to fill the vial to 10 ml to within 0.1 ml. This leads to
# errors for the grain density measurements on the order of 5 to 10 percent.

# Error propogation is done using gaussian error propogation

# As a check on our method, we can calculate the volume of the spherical beads accurately by measuring 
# the radius (2.475 mm) and calcluating the volume as 4/3*pi*2.475^3 = 6.35e-8 m^3. 
# We can also obtain a very accurate estimate of the average volume of the natural grains from 
# micro CT scans, which gives 3.64e-8.
# The volumes estimated using the below method are 6.15e-8 for the glass spheres and 3.53e-08 for the 
# natural grains. the ratio between our independent volume estimate and the volume estimate using
# the below method is 1.031 for the glass beads and 1.033 for the natural grains, suggesting that
# there is a systematic error where the grain volumes are estimated to be about 3% too small. This is 
# corroborated by the reported density of the glass beads of 2500 kg/m^3 which is lower than the estimate
# below by 5.5%. Therefore, we increase our estimate of the grain density below by 3% for all grain types.

# calculate density of grain assuming that the hole in several grain types is
# filled with water
def ro_p(grain_mass, grain_vol, hole_vol, grain):
    ro_w = 990.0
    print(grain, hole_vol)
    if hole_vol:
        return (grain_mass + ro_w*hole_vol) / (grain_vol + hole_vol)
    else:
        return grain_mass / grain_vol

# calculate error associated with density calculation using gaussian error
# propogation
def dro_p(grain_mass, dgrain_mass, grain_vol, dgrain_vol, hole_vol, dhole_vol):
    ro_w = 990.0
    if hole_vol:
        dpdm = 1/(grain_vol + hole_vol)
        dpdvolm = -(grain_mass + hole_vol*ro_w) / (grain_vol + hole_vol)**2
        dpdvolh = ro_w/(grain_vol + hole_vol) - (grain_mass + ro_w*hole_vol)/(grain_vol + hole_vol)**2
        return np.sqrt((dpdm*dgrain_mass)**2 + (dpdvolm*dgrain_vol)**2 + (dpdvolh*dhole_vol)**2)
    else:
        dpdm = 1/grain_vol
        dpdvolm = -grain_mass / grain_vol**2
        return np.sqrt((dpdm*dgrain_mass)**2 + (dpdvolm*dgrain_vol)**2)

# water volume and uncertainty on water volume in vial used for measurements
water_vol, dwater_vol = 1e-5, 1e-7 # m^3
dmass = 1e-6 # kg - error in mass measurements due to scale
rof, drof = 990.0, 5.0 # kg/m^3 - known density and uncertainty of water used
nu, g = 1e-6, 9.8 # kinematic viscosity of water and acceleration due to gravity

# make empty dictionaries to be filled with observations
empty_vial_mass, vial_with_grains_mass, vial_with_water_mass, vial_with_water_and_grains_mass  = {}, {}, {}, {}
grains, n_grains, ABC, hole_props = {}, {}, {}, {}

### Raw data
# Glass beads ('gb')
grains['gb'] = 'gb' # name grains
empty_vial_mass['gb'] = np.array([19.186, 19.179, 19.181]) / 1e3 # kg
vial_with_grains_mass['gb'] = np.array([22.427, 22.421, 22.417]) / 1e3 # kg
vial_with_water_mass['gb'] = np.array([29.069, 29.074, 29.122]) / 1e3 # kg
vial_with_water_and_grains_mass['gb'] = np.array([31.105, 31.124, 31.083]) / 1e3 # kg
n_grains['gb'] = np.array([20, 20, 20])
ABC['gb'] = None
hole_props['gb'] = None

# Faceted ellipsoids ('oc')
grains['oc'] = 'oc' # name grains
empty_vial_mass['oc'] = np.array([13.935, 13.928, 13.925]) / 1e3
vial_with_grains_mass['oc'] = np.array([16.390, 16.394, 16.368]) / 1e3
vial_with_water_mass['oc'] = np.array([23.911, 23.853, 23.869]) / 1e3
vial_with_water_and_grains_mass['oc'] = np.array([25.343, 25.339, 25.336]) / 1e3
n_grains['oc'] = np.array([15, 15, 15])
ABC['oc'] = None
hole_props['oc'] = {'hole_l': 0.005, 'dhole_l': 0.0001, 'hole_r': 5e-4, 'dhole_r': 1e-4}

# Natural river gravel ('ng')
grains['ng'] = 'ng' # name grains
empty_vial_mass['ng'] = np.array([19.207, 19.281, 19.186, 19.179]) / 1e3
vial_with_grains_mass['ng'] = np.array([25.010, 21.203, 22.214, 21.800]) / 1e3
vial_with_water_mass['ng'] = np.array([29.123, 29.058, 29.104, 29.127]) / 1e3
vial_with_water_and_grains_mass['ng'] = np.array([32.687, 30.268, 30.977, 30.732]) / 1e3
n_grains['ng'] = np.array([65, 21, 35, 29])
ABC['ng'] = None
hole_props['ng'] = None

# Cubic prisms ('ls')
grains['ls'] = 'ls' # name grains
empty_vial_mass['ls'] = np.array([13.926, 13.930, 13.928]) / 1e3
vial_with_grains_mass['ls'] = np.array([15.124, 15.188, 15.137]) / 1e3
vial_with_water_mass['ls'] = np.array([23.889, 23.884, 23.871]) / 1e3
vial_with_water_and_grains_mass['ls'] = np.array([24.581, 24.603, 24.579]) / 1e3
n_grains['ls'] = np.array([15, 15, 15])
ABC['ls'] = np.array([[3.5, 3.5, 4],[3.5, 4, 6],[3, 3, 4],[3.5, 3.5, 4],[3.5, 3.5, 5],[3.5, 3.5, 5],[3, 3, 3.5],[3.3, 3.3, 3.5],[3.5, 3.5, 4],[3.5, 3.5, 5]])
hole_props['ls'] = {'hole_l': ABC['ls'][:,2].mean() / 1e3, 'dhole_l': ABC['ls'][:,2].std()/np.sqrt(ABC['ls'].shape[0]) / 1e3, 'hole_r': 5e-4, 'dhole_r': 1e-4}

# Glass chips ('gg')
grains['gg'] = 'gg' # name grains
empty_vial_mass['gg'] = np.array([19.177, 19.179, 19.182]) / 1e3
vial_with_grains_mass['gg'] = np.array([23.872, 26.885, 25.152]) / 1e3
vial_with_water_mass['gg'] = np.array([29.030, 29.033, 29.083]) / 1e3
vial_with_water_and_grains_mass['gg'] = np.array([31.899, 33.720, 32.671]) / 1e3
n_grains['gg'] = np.array([20, 35, 25])
ABC['gg'] = np.vstack([np.array([[5.5, 7, 13.5],[3, 5, 6],[3, 7, 7],[3, 7, 9],[5.5, 6.5, 9],[5.5, 7, 11.5],[3.5, 4.5, 8],[3, 5.5, 7],[4, 6.5, 9],[3.5, 6.5, 8.5],[3.5, 6, 7.5],[3, 5, 7]]), np.array([[3, 5, 8], [5, 7, 13], [5, 6, 11], [4, 5, 12], [3, 4, 7], [2, 6, 12], [3, 7, 8], [4, 7, 10], [2, 5, 7], [5, 6, 7], [3.5, 6, 10], [2.5, 5, 6], [2, 6, 10], [4, 4, 10], [3, 5, 10], [3, 4, 5], [3, 4, 6], [3.5, 5, 7]])])
hole_props['gg'] = {'hole_l': ABC['gg'][:,1].mean() / 1e3, 'dhole_l': ABC['gg'][:,1].std()/np.sqrt(ABC['gg'].shape[0]) / 1e3, 'hole_r': 5e-4, 'dhole_r': 1e-4}


### Measuring the water density gives some idea of the accuracy and precision of
# the method we used
water_mass = {grain: vial_with_water_mass[grain] - empty_vial_mass[grain] for grain in grains}
dwater_mass = dmass
water_ro = {grain: water_mass[grain] / water_vol for grain in grains}
dwater_ro = {grain: np.sqrt((dwater_mass/water_mass[grain])**2 + (dwater_vol/water_vol)**2) * np.abs(water_ro[grain]) for grain in grains}
# print('Measurement error (in percent)', {grain: 100*np.abs(1 - water_ro[grain].mean() / rof) for grain in grains})
# print('Predicted error (in percent)', {grain: 100*(dwater_ro[grain].mean() / rof) for grain in grains})

# mass of water remaining in vial when grains are in vial and uncertainty
water_mass_exp = {grain: vial_with_water_and_grains_mass[grain] - vial_with_grains_mass[grain] for grain in grains}
dwater_mass_exp = dmass
# volume of water in vial when grains are also in vial and uncertainty
water_vol_exp = {grain: water_mass_exp[grain] / rof for grain in grains}
dwater_vol_exp = {grain: np.sqrt((dwater_mass_exp/water_mass_exp[grain])**2 + (drof/rof)**2) * np.abs(water_vol_exp[grain]) for grain in grains}

# grain volume (assuming total volume is 10 ml in vial) and uncertainty
# grain volume corrected by 3 percent to match manufacturer provided density of glass spheres.
grain_vols = {grain: (water_vol - water_vol_exp[grain]) / n_grains[grain] * 1.03 for grain in grains}
dgrain_vols = {grain: np.sqrt((dwater_vol)**2 + (dwater_vol_exp[grain])**2) / n_grains[grain] for grain in grains}
# grain mass in vial (using known mass of vial) and uncertainty
grain_masses = {grain: (vial_with_grains_mass[grain] - empty_vial_mass[grain]) / n_grains[grain] for grain in grains}
dgrain_masses = {grain: np.sqrt((dmass)**2 + (dmass)**2) / n_grains[grain] for grain in grains}
# grain density (using measured mass and volume of grain sample in vial) and uncertainty for each of three measurements
grain_ros = {grain: grain_masses[grain] / grain_vols[grain] for grain in grains}
dgrain_ros = {grain: np.sqrt((dmass/grain_masses[grain])**2 + (dgrain_vols[grain]/grain_vols[grain])**2) * np.abs(grain_ros[grain]) for grain in grains}
# grain hole properties
hole_vol = {grain: np.pi*(hole_props[grain]['hole_r']**2)*hole_props[grain]['hole_l'] if hole_props[grain] else None for grain in grains}
dhole_vol = {grain: np.sqrt((2.*np.pi*hole_props[grain]['hole_r']*hole_props[grain]['hole_l']*hole_props[grain]['dhole_r'])**2 + (np.pi*hole_props[grain]['hole_r']**2*hole_props[grain]['dhole_l'])**2) if hole_props[grain] else None for grain in grains}
# grain density accounting for water-filled hole in certain grain types and uncertainty for each of three measurements
total_grain_ros = {grain: np.array([ro_p(grain_masses[grain][i], grain_vols[grain][i], hole_vol[grain], grain) for i in range(3)]) for grain in grains}
total_grain_ros['ng'] = grain_ros['ng']
dtotal_grain_ros = {grain: np.array([dro_p(grain_masses[grain][i], dgrain_masses[grain][i], grain_vols[grain][i], dgrain_vols[grain][i], hole_vol[grain], dhole_vol[grain]) for i in range(3)]) for grain in grains}

### average across all three measurements of each grain type taking the average uncertainty as a measure of error for each sample of three
# average grain vol
grain_vol = {grain: grain_vols[grain].mean() for grain in grains}
dgrain_vol = {grain: dgrain_vols[grain].mean()/np.sqrt(grain_vols[grain].size) for grain in grains}
# average grain mass
grain_mass = {grain: grain_masses[grain].mean() for grain in grains}
dgrain_mass = {grain: dgrain_masses[grain].mean()/np.sqrt(grain_masses[grain].size) for grain in grains}
# average grain density
grain_ro = {grain: grain_ros[grain].mean() for grain in grains}
dgrain_ro = {grain: dgrain_ros[grain].mean()/np.sqrt(dgrain_ros[grain].size) for grain in grains}
# average grain density accounting for hole across three measurements for each grain type and uncertainty using standard error
total_grain_ro = {grain: total_grain_ros[grain].mean() for grain in grains}
dtotal_grain_ro = {grain: dtotal_grain_ros[grain].mean()/np.sqrt(dtotal_grain_ros[grain].size) for grain in grains}
# grain volume not including the hole
total_grain_vol = {grain: grain_vol[grain] + hole_vol[grain] if hole_vol[grain] else grain_vol[grain] for grain in grains}
dtotal_grain_vol = {grain: np.sqrt(dgrain_vol[grain]**2. + dhole_vol[grain]**2) if dhole_vol[grain] else dgrain_vol[grain] for grain in grains}

# save data to pandas dataframe and then to csv files.
df = pd.DataFrame({
                'water_mass': water_mass, 'dwater_mass': dwater_mass,
                'water_ro': water_ro, 'dwater_ro': dwater_ro,
                'water_mass_exp': water_mass_exp, 'dwater_mass_exp': dwater_mass_exp,
                'water_vol_exp': water_vol_exp, 'dwater_vol_exp': dwater_vol_exp,
                'grain_vols': grain_vols, 'dgrain_vols': dgrain_vols,
                'grain_masses': grain_masses, 'dgrain_masses': dgrain_masses,
                'grain_ros': grain_ros, 'dgrain_ros': dgrain_ros,
                'total_grain_ros': total_grain_ros, 'dtotal_grain_ros': dtotal_grain_ros,
                'grain_mass': grain_mass,'dgrain_mass': dgrain_mass,
                'grain_vol': grain_vol,'dgrain_vol': dgrain_vol,
                'hole_vol': hole_vol,'dhole_vol': dhole_vol,
                'grain_ro': grain_ro, 'dgrain_ro': dgrain_ro,
                'total_grain_ro': total_grain_ro,'dtotal_grain_ro': dtotal_grain_ro,
                'total_grain_vol': total_grain_vol,'dtotal_grain_vol': dtotal_grain_vol})

# df.to_csv('var_density_full.csv')

df2 = pd.DataFrame({
                'grain_mass': grain_mass,'dgrain_mass': dgrain_mass,
                'total_grain_ro': total_grain_ro, 'dtotal_grain_ro': dtotal_grain_ro,
                'grain_vol': grain_vol,'dgrain_vol': dgrain_vol})

df2.to_csv('1_var_density.csv')

# %%
